1. Field of the Invention
This invention relates to electrophysiological activity estimation methods which estimate equivalent current dipoles (i.e., ECD) based on the electromagnetic field distribution which is measured on a surface of a living body such as a scalp of a human. This invention is based on patent application No. Hei 9-19619 filed in Japan, the content of which is incorporated herein by reference.
2. Prior Art
Conventionally, the equivalent current dipole estimation method is known as the method to estimate active areas of a human brain on the basis of electromagnetic field distribution measured on a scalp of a human. By observing active areas of the human brain, it is possible to obtain information with regard to high-order functions of the brain as well as diseased parts of the brain.
Now, a general method to calculate an equivalent current dipole will be described with respect to estimation of active areas of the brain. When a nerve cell (or neuron) of the brain is excited by an impulse (or stimulus) given thereto from an external device, pulse-like current flows across a connection (or axon) connecting neurons together. Due to occurrence of the current, weak electromagnetic field is caused to occur around a scalp of the human. A source of the electromagnetic field is subjected to modeling using current elements called equivalent current dipoles (simply called "dipoles"). By estimating six parameters of the dipole, it is possible to estimate an active area of the human brain. Herein, the six parameters contain three parameters for designation of an area, two parameters for designation of a direction and one parameter for representation of magnitude (or intensity) with respect to the dipole.
Two methods are known as the method to estimate the parameters of the dipole. Herein, the six parameters of the dipole are unknown while a measured value of the electromagnetic field is represented by an equation as follows: EQU y=(y.sub.1,y.sub.2, . . . , y.sub.n).sup.T [Equation 1]
A theoretical value is calculated using the dipole model represented by an equation as follows: EQU .function.(.theta..sub.j)=(.function..sub.1 (.theta..sub.j), .function..sub.2 (.theta..sub.j), . . . , .function..sub.n (.theta..sub.j)).sup.T. [Equation 2]
So, a square residual is calculated as follows: ##EQU1## The aforementioned method estimates the parameter minimizing the square residual, which is represented by an equation as follows: EQU .theta..sub.j =(x.sub.j ',e.sub.j,q.sub.i).sup.T [Equation 4]
In the aforementioned equations, "n" denotes a number of measuring points; "m" denotes a number of dipoles; "y.sub.i " denotes a measured value of the electromagnetic field measured at a measuring point i; ".function..sub.i " denotes a theoretical value of the electromagnetic field at the measuring point i. In addition, a vector corresponding to combination of the measured value and theoretical value of the electromagnetic field is represented as follows: EQU y,.function.
Incidentally, an expression of (. . .).sup.T represents the transposition. Further, a coordinate of the measuring point i is represented by EQU x.sub.i
while a dipole parameter at the measuring point j is represented by an equation as follows: EQU .theta..sub.j =(x.sub.j ',e.sub.j,q.sub.i).sup.T [Equation 5]
In addition, the position, direction and magnitude of the dipole at the measuring point j are respectively represented by symbols as follows: EQU x.sub.j ',e.sub.j,q.sub.j
As the method to estimate the foregoing dipole parameters, it is possible to use the nonlinear optimization algorithms such as the Levenberg-Marquardt method and simplex method. Details of the nonlinear optimization algorithms are shown by "paper 1" which is a book entitled "Nonlinear Programming" which is written by Mr. Hiroshi Yamashita and Mr. Hiroshi Konno and is provided by Science Technology Association of Japan on 1987, for example. Further, the method using the pseudoinverse of the matrix (or Moore-Penrose inverse) is known as an example of the other methods for the current dipole estimation, especially as an example of the estimation method for "distributed" active areas of the living body. According to this method, the theoretical equation of the electromagnetic field is represented as follows: ##EQU2## So, it is possible to use the characteristic of the above equation that elements regarding the magnitude of the dipole are coordinated by the linear function. Concretely speaking, the shape of the brain is approximated by the set of the polyhedrons while the position and direction of the dipole are fixedly located on the polyhedron, so that the magnitude q.sub.i of the dipole is calculated as an unknown parameter. In this case, the function of "F" is represented by the matrix, so the aforementioned equation 6 is rewritten in the matrix form as follows: EQU .function.=Fq [Equation 7]
For example, the electroencephalograph (i.e., EEG which represents a brain wave measuring instrument) is used to measure the electric field distribution (or electric potential distribution), based on which the dipole is estimated. In this case, calculation for F is made with regard to elements i, j in accordance with an equation as follows: ##EQU3## In the above equation, "Y" denotes the spherical harmonics; symbols r,.theta.,.phi. show the position of the measuring point in polar coordinates while symbols r',.theta.',.phi.' show the position of the dipole in polar coordinates. A symbol .gradient.' denotes a differential operator while a mathematical expression of R(r,r') corresponds to the function regarding r and r', which is determined by the boundary condition. Using a SQUID (an abbreviation for "Superconducting Quantum Interference Device"), it is possible to measure the magnetic field distribution, based on which the dipole is estimated in accordance with an equation as follows: ##EQU4## where .mu..sub.0 is a magnetic permeability. In addition, n.sub.i denotes a normal vector with respect to a surface of a head model at a measuring point i. A square error is represented by an equation as follows: EQU E=.parallel.y-Fq.parallel..sup.2 [Equation 10]
Herein, a vector q which minimizes the above square error is represented, using a pseudoinverse F.sup.+ of F, by an equation as follows: EQU q=F.sup.+ y [Equation 11]
The dipole estimation method using the pseudoinverse is disclosed by "paper 2" corresponding to Japanese Patent Laid-Open Publication No. 6-343613 as well as "paper 3" corresponding to Japanese Patent Laid-Open Publication No. 6-343614, for example.
Now, FIG. 6 is a flowchart showing an example of processing which the conventional system executes in accordance with the conventional method using the pseudoinverse. First, a group of grid points are arranged uniformly with respect to a diagnosing area, i.e., an area on which active areas are estimated in step 51. In step 52, a current source is calculated with regard to each of the grid points in accordance with the minimum-norm least-squares method. In step 53, a group of grid points are moved in proximity to a grid point on which a current source having a large value exists. Then, a minimum space between grid points is selected from among spaces between the grid points. In step 54, a decision is made as to whether the minimum space is less than a prescribed value or not. If the minimum space is less than the prescribed value, the system completes the proceeding of FIG. 6 because the convergence is obtained. If not, the system proceeds back to step 51.
According to the aforementioned method, the system performs estimation while changing the distance between the grid points in order to improve a precision of estimation. However, the conventional system does not perform estimation with respect to a number of dipoles. In addition, the conventional system is not designed to cope with problems due to the noise.
Next, a description will be given with respect to problems which the conventional dipole estimation method suffers from.
Normally, the method using the nonlinear optimization algorithm uses one to three dipoles to perform estimation of an active area. However, it is impossible to know a number of the active areas in the brain in advance. For this reason, there are provided dipole models whose number of dipole(s) is one, two and three respectively, for example. Using the dipole models, the system performs the estimation. So, the system employs the best result of the estimation. In this case, a square error E can be used as an element of evaluation criterion for the quality of estimation. However, there is a tendency that the square error becomes smaller as a number of dipoles is increased more. For this reason, there is a problem that a number of dipoles cannot be detected with accuracy. FIG. 7 shows an example of variations of the square error of the estimation which is made with respect to the magnetic field distribution produced by two dipoles while changing a number of model dipole(s) from one to eleven. In FIG. 7, the vertical axis corresponds to the logarithmic notation of the mean square error (or mean square residual, i.e., MSR) while the horizontal axis represents a number of model dipoles. According to the graph of FIG. 7, it is observed that the square error is decreased as a number of model dipoles is increased. So, if the square error is used as the evaluation criterion, although a true number of dipoles is two, the best result of the estimation is obtained with respect to the case of eleven model dipoles. This is a wrong result. In addition, the method using the nonlinear optimization algorithm requires a long time to perform the estimation one time. Further, it is necessary to perform the estimation several times while changing the number of the dipoles. Namely, the conventional method requires a large load to calculations as well as a long time for the estimation.
Another conventional method using the pseudoinverse arranges several hundreds to 10,000 dipoles with regard to each element of the polyhedrons approximating the shape of the brain, wherein magnitude of each dipole is calculated. According to this method, the estimation can be made using multiplication of matrix only. So, this method requires a short time for the estimation. Even if a neural activity lies in an extended area, the method is capable of providing information about the size of the area. The method using the pseudoinverse is designed to use the linear model. For this reason, this method has a drawback corresponding to weakness against the noise.
Now, a description will be given with respect to reasons why the method using the pseudoinverse has such a drawback.
Using the singular value decomposition (i.e., SVD), the pseudoinverse F.sup.+ is written by an equation as follows: ##EQU5## where u.sub.i,v.sub.i represent column vectors i for matrices U, V; .lambda.i represents a singular value of the matrix F, so sum is performed with respect to singular values which are not zero. If measurement data do not contain noise, correct dipole distribution can be obtained by multiplication of the pseudoinverse and measurement data in accordance with an equation as follows: ##EQU6## where q corresponds to a solution vector having a minimum norm which minimizes the square error. However, in the case where the measurement data contain noise z, the equation 13 should be rewritten by using an equation as follows: EQU y=y.sub.0 +z [Equation 14]
where y.sub.0 corresponds to measurement data containing no noise. So, using the equation 14, the equation 13 is rewritten as follows: ##EQU7## In the above equation, a second term of the right side shows an extremely large value with respect to a small singular value .lambda.i. So, the divergence occurs on an estimated value of the magnitude of the dipole.
So, the problems of the conventional methods can be summarized as follows:
In the case of the method using the nonlinear optimization algorithm where the square error is used as the evaluation criterion for the estimation, it is impossible to estimate a number of dipoles. In the case of the method using the pseudoinverse for estimation of the distributed active area, the divergence occurs on estimation of the magnitude of the dipole due to the effect of the noise.